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I. INTRODUCTION 


For a structural analysis of any non-trivial elasto-static 
problem by the finite element method, a large amount of data must 
be provided to the computer system which performs the calculations. 
The required information includes the arrangement of nodes within 
each element, the coordinates of all nodes, and the applied loads. 
After the numerical values of this data have been determined, which 
is a large problem in itself, the information must be supplied to 
the computer system. This is normally accomplished by punching the 
appropriate data on computer cards. The typing of these cards is 
also a difficult and time-consuming iob. More important, however, 
is the ever-present vossitiiity of numan errors during the performance 
of these tasks. If an error is not detected prior to submitting 
the data for eeeeeeee expensive computer time would be wasted. 

With the above considerations in mina, it is readily apparent 
that a computer system, which would generate most of the data required 
by the primary computer system for a complete structural analysis, 
would be very useful. In addition, this mesh generator should 
require as little input data as possible. 

References 5, 8, and 10 contain reports of efforts which have 
been directed at coding computer systems which will aid the analyst 
in the interpretation of output data. A few computer svstems, which 
are reported in Refs. 1, 2, and 9, have been coded to automatically 


generate input data. However, these input data generation systems 





are either for types of finite elements other than isoparametric 
elements or for two-dimensional elements. 

Two input data generator systems for chteosdinencieeln iso- 
parametric finite elements will be presented in this thesis. One 
System is for quadratic elements and the second system is for cubic 


elements. 





Ti, ILSOPARAMETRIC FINITE ELEMEMiS 


The isoparametric concept will be discussed in this chapter in 
sufficient detail to establish only the procedures for constructing 
isoparametric elements and the validity of relationships used to 
obtain the consistent gravity and pressure loads. The readers are 
referred co Refs. 4 and 7 for further details of the isoparametric 


concept. 


A. DISPLACEMENT FUNCTIONS 
The displacements of nodes in any one element are obtained by 
uSing shape functions which are defined in a system of dimensionless 
coordinates (§,7.C) that are unique for that element. The coordinates 
range from minus one to plus one. 
= N 
u(x, y, 2) = YN.(§,1,0) u, 
TACs nue aN (5), 516) V. 
1 i: 1 


aly. -¥ selZ) o> % NS ,1,6) we 


(1) 


where N. 51,6) are the shape functions and Uns Vis W, are the dis- 


placement components of node "i" in a global reference system (x, y, 2), 


where 'i'' assumes values from 1 to the total number of nodes in the 


element. This global reference system is then related to the dimen- 
Sionless coordinates by the same shape functions as follows: 
x(§,1,0) = ON. G,1,0) x; 
ot 
y(§,1,6) = DN, (8.1,0) yy 
a: 


2(6,,6) = DN, (81,0) 2 


(2) 
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where xX, Y,» 2, are the coordinates of node "i'' in the global 
reference system [4]. 
The shape functions referred to abov2 have been obtained from 


(7]. They are Listed in section 2.C, 


B. TRANSFORMATIONS 

In order to obtain the consistent gravity and pressure loads, 
the volumetric increment, dxdydz, must be found in terms of the 
dimensionless coordinates. That operation requires the Jacobian 
matrix [4]. 


ax/d§ dy/a§ az/ag 


{s] = ax/dT dy/dt] az/an (3) 
ax/dC ay/ot az/ac 
The element of volume is then transformed as follows: 
dv = det{J|d& dt dc ; (4) 


C. SHAPE FUNCTIONS 

In this listing of the shape functions, attention is directed 
to the tri-quadratic and tri-cubic elements, for which the mesh 
generator programs were written. | 

1. Quadratic Elements 

see Figure 1 for identification of the nodes. 

Maemer nodes: 1, 3, 5, 7, 13, 15, 17, and 19. 

ML, CUE) GUes Ie NCL rip yGcr 6 3G aries Gaal) Gs) 
where on = 5,5 and 24) rife ia ee sant lan eyo rs No and Coe 
Mid-side nodes: 2, 6, 14, and 18. 


_ 2 
N, = (1/4) (1 = 6) + 1) + 6) (6) 


ll 





Mideside nodes¢ms4. 6. .'6- anca.20, 
9 
N, = (U/L = WA +80 +6) (7) 
Mid-side nodes: 9, 10, 11, and 12. | 


2 
Nees G4) Cla a S)) Cae ie (8) 


2, Cubic Elements 
See Figure 2 for identification of nodes, 

Morner nodes: ~ 15°45 75° 10, 21,24, 27 mane 20, 

N, = (1/64)(L +8) +11 + C096" + 1° +67) ~ 19) (9) 
Mid-side nodes: 2, 3, 8, 9, 22, 23, 28, and 29. 

N, = (9/64) - BA 49H VA H+TIAFC) (10) 
where § = 7 1/3, 1, = 21, and ¢, = 91. 
Mid-side nodes: 5, 6, 11, 12, 25, 26, 31, and 32. 

N, = (9/64)(L = T°) + MA +8) +6) ab 
where N. = B38 55 =o i. Gi =f Ibe 
Mid-side nodes: 13, 14, 15, 16, 17, 18, 19, and 20. 

N, = (9/4 (L- CVA + DA +E) +1) (12) 


O 
where C. == lye S. = 71, and ie me 


Note that the convention for numbering of nodes within any 
element is to begin at the location where § = ]| = € = 1 and to 
consecutively number the nodes, proceeding in a counter-clockwise 


direction around the C-axis (see Figures 1 and 2). 


12 





Figure L. Quadratic Element 
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Figure 2. Cubic Element 
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IIL. DISCUSSION OF COMPUTER PROGRAMS 


The computer programs presented in this thesis were written to 
Support a computer system called TRISOP. TRISOP was coded by Professor 
G. Cantin of the U. S. Naval Postgraduate School, Monterey, California. 
TRISOP performs a structural analysis of three-dimensional elasto- 
Static problems using isoparametric finite elements. Information 
required by TRISOP for any problem includes general mesh parameters, 
title cards, format statements, and material properties; Young's 
modulus of elasticity, Poisson's ratio, and the coefficient of thermal 
expansion. Fifteen cards are needed to supply this data to TRISOP. 

The bulk of the data required by TRISOP is the element connectivity, 
structure geometry, loading, and boundary conditions. The boundary 
conditions will vary with each problem; therefore, no effort has 
been made to automate this process. 

The mesh generating computer programs discussed in this chapter 
compute the remainder of the data required by TRISOP: element 
connectivity, coordinates of nodal points, consistent gravity load 
vector, and consistent pressure load matrix. Additional information 
which rs helptul ime preparine an input deck of «cards for) [RfSOPeis 
also computed. One of these programs generates data for quadratic 
elements and the second program generates data for cubic elements. 


These programs will be called QUAMEG and CUMEG. 


A, ELEMENT CONNECTIVITY 
A structure of any shape can be visualized as a cubical piece 


of plastic material wnich can be stretched, compressed, and otherwise 
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appropriately deformed in such a way that it will assume the geometry 
of the structure. The element connectivity is determined when the 
structure is in this initial abstract cubical shape. The dimensionless 
coordinates of this abstract structure are & , tee Car 

Numbering of nodes is begun in a face which will have the 
fewest nodes, In that face, numbering proceeds sequentially along 
the edge having the smallest number of nodes as shown in Figure 3. 
This convention reduces the size of the half-bandwidth of the system, 
which directly affects the space and time required by the computer 
to solve the problem. The coordinate system is positioned such 
that node one is at ate = 0, ie = 1, and pe = 1. An example is 
discussed in Appendix A. 

For this mesh, the connectivities of the elements follow simple 


but tedious rules that were coded directly. 


-~ 


B. SOLUTION TIME AND SPACE REQUIREMENTS 
The time required to solve the system of equations in the problem 

and the space required by the computer to store all data are important 
information. They must be determined before the problem is submitted 
for processing. The caiculation of this information requires no 
additional input data and it is therefore included in QUAMEG and CUMEG 
to reduce the effort required of the programmer. 

The time required to solve the equations is computed using 
an empirical formula which was constructed by Cantin [4]. This 
calculation has been included in the coding of QUAMEG and CUMEG. 

The space required to store the various bits of information 


for a problem depends upon the total number of records and the 
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length of each record. These will vary with di fferentee pesuan 
computers. Simple formulas which generate the correct results for 
the IBM 360/67 computer used at the U. S. Naval Postgraduate School 


are employed in the coding of QUAMEG and CUMEG, 


C. NODES ON BOUNDARY SEGMENTS 
Before the coordinates of nodes which are within the boundaries 
of a structure can be computed, the coordinates of all corner nodes 
on the edges of the bounding surfaces must be specified. This problem 
is simplified by the method discussed below. — 

The term "boundary segment'’ must first be defined for the purpose 
of this discussion. A segment nee satisfy a few simple conditions 
which depend upon the orientation of the proposed segment in the 
Structure and upon the ctype of coordinate system uSed to descrive 
the geometry of the structure, QUAMEG and CUMEG are coded to allow 
the use of either rectangular or cylindrical coordinate systems. A 
combination of the two coordinate systems in any one problem is not 
permitted. 

1. Rectangular Coordinates 

When rectangular coordinates are used, a segment is any 
Straight Line between two corner nodes of the mesh. The element 
divisions along this line must be uniform. There are no other restric- 
tions imposed upon a segment when rectangular coordinates are used. 

2. Cylindrical Coordinates 

When cylindrical coordinates are used, the radius and reference 
angle are taken in planes parallel to the x-y plane. In this case 


a segment can be an arc of circle centered on z and ina plane 
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parallel to plane x-y, or a straight line in a plane containing 
the z-axis. Again, element divisions along these lines must be 
uniform. 

In order to determine the coordinates of nodes along a boundary 
segment, the coordinates of the corner nodes at the ends of the 
segment must be provided in the input data. In addition, the location 
of the segment within the mesh must be specified. The method of 
specifying the location of a segment is described in sections 4.B.3 
and 4.B.4, 

In Figure 4, the vector AB is easily obtained from the coordinates 
of nodes A and B. After calculating the values of increments of 
distances along x, y, and z, the intermediate nodes are easily 
obtained. 

If the bounding surfaces of a structure are of an irregular 
configuration such that portions of the edges of that surface cannot 
be specified by segments, the coordinates of all nodes along this 
irregular portion of the bounding surface must be provided in the 


input data. 


D. CORNER NODES 

Once the coordinates of the nodes on all of the boundary segments 
are established, the Poord nates of corner nodes on the bounding 
surfaces and on interior slices are calculated using a method described 
below. The bounding surfaces are considered in the following order 
(Bigure 3): face 1 (C’=1), face 2 (¢ =0), face 3 (E‘=0), bacers (E'=1), 
mace 5 (N’=1), face 6 (N=0) . Slices which are each of constant a 
are then considered until the coordinates of all corner nodes are 


calculated. 
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Figure 4. Boundary Segment 





The method used for calculating the coordinates of corner nodes 
within the edges of a bounding surface or slice can best be understood 
with the example illustrated in Figure 5. In this example, the 
surface ABCD is the real surface for which a mesh is desired. The 
mesh is first constructed in the a8 plane as shown. 

At this time, the coordinates of all the corner nodes on the 
boundary are known. If the coordinates of the interior nodes are 


required to satisfy the Laplacian equation: 


wae = 0 
v (y) =~ 0) 
V7*(z) = 0 Gis) 


where Th = 3° fac’ + ORS 

a set of acceptable interior nodes will be obtained. The dummy 
variables (@,8) are replaced by appropriate combinations of a a, 
and e depending upon the location of the face being considered. 
For this computation, the standard Laplacian molecule of finite 
difference theory is used and the final coordinates are obtained 

by relaxation. Since the coordinates of the nodes within the edges 
of a face are arbitrary, the relaxation process can be stopped a 
any time, depending upon the degree of approximation desired. In 


this code, sixty iterations are used. 


E. MID-SIDE NODES 

After the coordinates of all corner nodes have been established, 
the coordinates of the mid-side nodes are computed by averaging 
the coordinates of adjacent corner nodes. Although the mid-side 
nodes need not be Located exactly mid-way between the corner nodes, 


it is the simplest means of determining their coordinates, 


Za 








face 


Figure 9- Bounding Sur 
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F, CONSISTENT GRAVITY LOAD VECTOR 
The z-axis has been Shosen as the tine ot action of the gravity 
load. The work done by the force due to gravitational acceleration 


on a body is 


U = \ pgwdV (14) 


where U is the work done; Pp is the mass density of the body; & is 
gravitational Fecelerationemws 15 5c displacement of the body in 
the z-direction; and dV is the element of volume. 

The work done by concentrated forces acting at each node in 
the structure mesh is, in matrix notation: 

U = dug? (15) 
where (V) is the consistent gravity load vector and (w,) are the 
respective node displacements. 

From Chapter 2, recall the following two expressions: 

wey.z) = ENE. YMG 

av = det[ J] d§d dc. 

When these expressions are substituted into equation 14 and the 
result identified to equation 15, the consistent gravity load vector 


is obtained. 


lealel 
wy = pe \ VF aetlskn, ezanes (16) 
“Teil 


In this form, equation 16 ig ideally suited for solution by 
Guassian cubature. In QUAMEG, four Gauss points are used in each 
of the three directions, for a total of 64 points within each element. 
Five Gauss points are used in CUMEG. The weight factors Bheval (S@versc 
dinate values are obtained from Ref. ll. The Jacobian matrix 1s 


determined by using three function statements and the shape functions 
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are determined with two function statements in Subroutine GRAP 
(Appendix ©). These function statements are permuted to obtain the 
partial derivatives of the shape functions and the values of the 


shape functions for all nodes in the mesh. 


G. CONSISTENT PRESSURE LOAD MATRIX 

The consistent pressure load matrix for a uniformly distributed 
pressure acting on any face of an element is determined in a manner 
similar to that used for the gravity load. 


The work done by the uniformly distributed pressure is: 


u = \\(-pi da): @) (17) 


where n is a unit normal vector (positive outward) and U is the 
displacement vector at the surface (see Figure 6). The displacement 
vector u can also be written as: 


u(x,y,z)1 + v(x,y,z)4 + wx,y,z)k (18) 


c | 
N 


The vector ndA is found from differential geometry: 
ndA = (dr/dv)X(dr/a8) dads (19) 
where Y = xi + val +zk isa position vector of a point on the pres- 
surized surface and a and $8 are dummy variables, standing for either 
oer, or C. 

The terms dr/odd and dr/d8 are simply two rows-of the Jacobian 
matrix evaluated on the surface. After making appropriate substi- 
tutions and performing a few matrix manipulations, the expression 


meauces to” 


lL 1 
U = ( \ (£(v,8)) sth dads (20) 


-l -1 
where oe is the vector of nodal coordinates for the eight nodes 


on the face and £(@,8) is a row vector of functions of @ and 8B. 
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The integration is performed by Gaussian quadrature and the 
resultant vector is exhibited in convenient matrix form. Five Gauss 
points are used in both QUAMEG and CUMEG. 

The Jacobian matrix and shape functions are determined in 
Subroutine GRAP by using the same function statements as are used 
in the calculation of the consistent gravity load vector. The face 
number of the element on which the pressure is acting must be included 
in the input data. With this information available, a "computed 
go to'’ statement is used in the code to determine which variables 
(§, Tl, €) will replace @ and 8. Referring to Figure 6, the faces 


are located as described below. 


Face Location @ 8B Face Location ao 8B 
1 €=H 71 h 2B Si iy 
2 Uae Gs 5 evel S@ 4G 
3 C= +1 al 6 C =-1 My 


Wee LOT OF STRUCTURE MESH 

Once the coordinates of all nodes in the mesh have been established, 
the mesh must be inspected to determine if it is acceptable. Without 
a plot of the coordinates, this task is nearly impossible. Plotting 
the mesh by hand is prohibitively Seer Subroutine GRID 
has been included in QUAMEG and CUMEG to generate an off-line printer 
plot of the mesh. The parameters required for IBM 360 source library 
program DRAW are determined in GRID. DRAW is called to perform the 
actual plotting of the mesh. GRID generates a two-dimensional 
plot of the mesh on slices at constant values of Ge For the example 
shown in Figure 3, faces 1 and 2 and five slices at constant values 


of C° would be plotted. 
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The user can specify that the scaling be performed automatically 
by DRAW. However, that choice could result in the x-scale being not 
equal to the y-scale. The plot of the mesh would be distorted in 
that case. The scaling is done in GRID to inSure that the scales 
will be equal and that the smallest possible scale is used. 

The plotting space available with DRAW is fifteen inches in 
the y-direction and nine inches in the x-direction. Recall that 
the convention for numbering the nodes usuelly results in having 
the narrowest side of the structure along the y-axis (section 3.A). 
Therefore, GRID interchanges the x= and y-coordinates of the structure, 
for plotting purposes only, to take advantage of the larger plotting 
space in the y-direction. This further reduces the scale of the 


I 


plot. 
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IV. INPUT DATA PREPARATION 


This chapter is written as a self-contained unit giving specific 
information necessary for the preparation of input data cards. The 
chapter may be removed and used as a guideline for input data preparation 


after a thorough knowledge of the computer programs has been obtained. 


A. VARIABLES 
The variables and formulas defined below are included in the 
code and should be used as an aid for determining the arrangement 


of nodes within the mesh when preparing input data. 


Waagke) be) Meaning 
i Wil Total number of elements. 
NEL = (NX) (NY) (NZ), where NX, NY and NZ are the number of elements 
along the Zee and ies directions (21) 
Ze NPS Total number of nodes on any face at constant 


if ‘ C ° e 
G , where a face is a slice which intersects 
both corner nodes and mid-side nodes. 


NPS = 3(NX) (NY) + 2(NK + NY) +1 (QUAMEG ) (22) 
NPSe= S(NX ONY) + SCN + NY ee 1 (CUMEG) - §425)) 
oe NPM Total number of mid-side nodes between any two 
faces. 
NPM = (NX + 1) (NY + 1) - (QUAMEG ) (24) 
NPM = 2(NX + 1)(NY + 1) (CUMEG) (25) 
2 NPL Total number of nodes in one layer, where a 
layer consists of one face plus the slice 
intersecting mid-side nodes between that face 
and an adjacent face. 
NPL = NPM + NPS , (26) 
ae NUMNP Total number of nodes in the mesh. 
NUMNP = (NZ) (NPL) + NPS (27) 


28 





3. DATA CARDS 
i, Element Field Data (ii L072) 


Number of cards per problem - one 


Columns Variable Meaning 

i= 10 NX Number of elements im ocne € =direction, 
ee) 20 NY Number of elements gh So N\'-direction; 
21 - 30 NZ Number of elements *n the © -direction; 
31 - 40 ND X Number of segments iinet We Ee’ direction; 
41 - 50 NDY Number of segments Ime ewe 1\"-direction; 
51 - 60 NDZ Number of segments in the ¢’-direction; 
61 - 70 NPDP Total number of nodes for which pelnnsy  (eveloie= 


dinates must be specified to define the 
location of segment end points or other 
boundary nodes; 


71 - 75 KORD Specifies whether rectangular oF cy Lindri- 
cal coordinates 2Te used in the problem. 
A positive integer indicates cyiindricse! 
coordinates and a blank or zero specifies 
rectangular coordinates ; 


76 - 80 NCARD A positive integer indicates that punched 
cards are desired as output. A blank or 
zero indicates that cards are not desired. 


2. Problem Identification (10A8) 


Problen ica 


Number of cards per problem ~- one 


Columns Variable Meaning 
1 - 80 TITLE Alphameric information [to identify the 
problem. 


3. Nodal Point Coordinates (1110, 3F10.0) 


Nodal_toint 


Number of cards per problem - one for each NPDP point 


Columns Variable Meaning 
1 - 10 i Nodal point number 5 
11 - 20 COORD(1,1) x-coordinate (radius) ; 





Columns Variable Meaning 
21 - 30 COORD (1,2) y-coordinate (angle); 
31 - 40 COORD Cis3) z-coordinate. 

4. Segment Divisions (5110) 


Number of cards per problem - one for each segment 


Columns Variable Meaning 

1 - 10 iL Lowest node number in the segment ; 

11 - 20 J Highest node number in the segment; 
m= 30 N1 Segments in the Ze and (‘-directions: 


Number of elements above the segment 

(in the positive 1\"-direction); Segments 
in the 1) -direction: Number of elements 
within the segment; 


31 - 40 M1 Segments in the ee and C’-directions: 
Number of elements to the left of the 
segment (in the negative 6 ‘-direction); 
Segments in the ~~ d treeeren: He trlegers 
of elements within the segment; 


aie= 50 Ll Segments in the €’- and \’-directions: 
Number of elements before the segment 
(in the positive C‘-direction); Segments 
in the C'-direction: Number of elements 
within the segment. 


PeeboadiData V2 il0, ZE10, 0) 


Number of cards per problem - one 


Columns Variable Meaning 
1 - 10 MAP A positive integer indicates that a 


plot of the structure mesh 1s desired. 
A blank or zero indicates that a plot 
is not desired; 


11 - 20 NEFP Number of pressurized element faces; 
Zi= 30 SGZ Specific weight ; 
31 - 40 UDP Uniformly distributed pressure. 
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6. Pressure Loading Data (2110) 


Number of cards per problem - one £or each NEE 


Columns Variable Meaning 
1 - 10 NELP Element number on which a pressure exists; 
11 - 20 NFACE Face on which the pressure is applied. 


(See section 3.G and Figure 6) 
7. Plot Title Cards (6A8) 


Number of cards per problem - two per ey Ligtie (Oe a= Uh plots) 


Columns Variable Meaning 
1 - 48 EEE Alphameric information to identify the 
plot; 
1 - 48 ITITLE Alphameric information to POentr ey) eme 
programmer and his computer center box 
number. 


8. Stop Trap 


The stop trap employed in the programs ailows tae execution 
of as many separate problems as may be desired with one submission 
of the program. The program will continue processing problems Miia 
it is informed that NX is less than zero. Therefore, the last card 
of a complete input data deck should have a negative integer in 


columns 1 - 10. 
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Vv. CONCLUSIONS 


The addition of QUAMEG and CUMEG to the tools available in the 
finite element method of analysis opens new doors for the analyst 
who uses iscparametric finite elements. Previously, the amount 
of input data required to analyze a structure of appreciable size 
was too overwhelming for the analyst to undertake the problem. 
QUAMEG and CUMEG make it practical for the analyst to solve such 
a problem. In addition, he can be certain, after an acceptable 
mesh is generated, eae the output data generated by QUAMEG and 
CUMEG is free of error. 

For the example discussed in Appendix A, the number of input 
cards required by QUAMEG versus the number of output cards generateu 
by QUAMEG gor direct input to TRISOP is about six per cent. For 
the 1 x 8 x 8 mesh of the pinched cylinder problem discussed in 
Appendix B, this ratio is about 4 per cent. 

With a gravity load and pressure load acting on one cubic 
element, it was found that QUAMEG and CUMEG produced consistent 
load vectors which were accurate to at least eight decimal places. 

QUAMEG and CUMEG Bees de rabies reduce the effort required of 
the analyst to solve a problem by the finite element method, generate 
acceptable values for nodal point coordinates and compute accurate 


consistent gravity and pressure load vectors. 
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VI. RECOMMENDATIONS 


A. LOADING CONDITIONS 

The capabilities of QUAMEG and CUMEG should be increased by 
adding subroutines which will generate the consistent load vectors 
for a thermal gradient in the structure and for an initial displace- 


ment condition at one or more nodes. 


peeeeLOT OF STRUCTURE MESH 

It would be preferred to obtain a more complete representation 
of the structure mesh than is now provided. Two courses of action 
are recommended: 

lL. &xpand Subroutine GRID in such a way that plots of faces 
in the eee. and/or 2 5c planes are also generated. 

2. Replace GRID with a subroutine which would generate a 
three-dimensional plot of the structure mesh. The feasibility 
of using the IBM 360 source library program PLOT PACKAGE, which 
is mentioned in Ref. 12, should be investigated if this method 
is attempted. 

This author feels that the results of recommendation 2 would 


be preferred to 1, although the coding may be more difficult. 
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APPENDIX A 
GENERATING A MESH WITH QUAMEG 
This appendix feeeetbes the preparation of input cards and the 
expected results for a sample problem using QUAMEG. This problem 
was chosen as an example to demonstrate the capabilities of QUAMEG 
and to more fully explain the procedures described ia Chapter 4. 

The structure considered is a cylindrical solid with a star 
shaped void in the center. This void linearly converges to a point. 
Figure 7 shows a 45° sector of the ee ee It is normally prefer- 
able to work with angles in the first quadrant. Therefore, the 
sector from zero.to forty-five degrees will be considered. 

The first step in ihe problem is ta decide upon a mesh size. 

In this case, an 8 x 5 x 4 mesh was chosen. The amount of control 

to be exercised by the programmer to inSure that the structure 

will be modeled accurately can be determined now. The coordinates 

of all nodes along the edges of the bounding surfaces of the structure 
must be available to the system before the remaining nodal coor- 
dinates can be determined. Since there are no irregularities on 

the boundaries between the front and back faces, the coordinates 

of the nodes on the edges between these two faces can be computed 

by providing the necessary segment data. Therefore, only the nodes 

on the front and back faces need be considered. 

At this point, two rectangles of constant ie are drawn with 
the proposed mesh superimposed upon them (ee Figures 8 and 9). 


These figures will serve as an aid for numbering the nodes on the 
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actual structure and for determining the locations of segments. 
After the nodes on face one are numbered, equation 22 is solved to 
insure that an error has not been made in the numbering. The result 
of equation 22 is subtracted from that of equation 27 to determine 
the number of the first node on the back face. In this way, the 
programmer need not go through the time-consuming process of numbering 
all nodes in the mesh. The nodes on the back face are numbered and 
checked as before. Note that the 8 x 9 mesh is in the eae plane. 
Orienting the structure so that the 5 x 4 face of the mesh were in 
the eel, plane would reduce the half-bandwidth of the system. 
However, cylindrical coordinates must be used in this problem and 
it would not be possible to model the structure accurately if ltmwere 
oriented in’that manner. Finally, a picture is drawn of the front 
and back faces of the structure with the desired mesh super iwposce 
and the nodes numbered (Figures 10 and jie 

The identification of segments will now be determined. The 
segments will be considered in the same order as they are considered 
in the program. Referring to Figures 8 and 9, it is seen that the 
following lines are parallel to the Neaxis: 1 - 11, 137 - 147, 
805 - 815, 941 - 951. Now referring to Bicures 0 andsn ll.) aes 
seen that all of these lines are of constant angle with a uniform 
distribution of elements within the lines. Therefore, these lines 
will be treated as segments and the coordinates of the nodes listed 
above will be included in the input data. 

By this method of specifying lines, the author does not intend 
to imply that the nodes along the line are numbered consecutively, 


starting with the lower number. For example, the nodes along ine 
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60 - 111 in Figure 8 are: 60, 67, 77, 84, 94, 101, and I11. This 
method of specifying lines merely lists the numbers of nodes at the 
ends of a Pine: 

In order to control the departure angles of lines 79 - //, 

96 - 94, and 113 - 111, it is necessary to force line 52 - 62 nearer 
the x-axis than would result if an even distribution of nodes along 
line 1 - 137 were allowed. Therefore, an appropriate location for 
node 52 is chosen, the coordinates of nodes 52 and 62 are included 
in the input data, and line 52 - 62 is treated as a segment. Note 
that this line is of neither constant angle nor constant radius. 
However, Since the line is not on the edge of a bounding surface, 
the fact that it will be curved will not affect the solution adversely. 
Similarly, the coordinates of nodes 856 and 866 will be included 

in the input data and that line will be a segment. 

The same rules govern the selection of segments in the €'-direc- 
tion. Referring to Figures 8 and 9, it 1s seen that the coordinates 
of nodes along lines 1 - 137, ll - 147, 805 - 941, and 815 - 951 
must be made available to the computer. In Figures 10 and il, line 
1 - 137 is a line of constant radius; however, the element divisions 
along this line are not uniform. Two segments must be specified: 
Memes 1 ~- 52 and 52 = 137. The coordinates of the end nodes on these 
segments have been specified already and nee rot be repeated. Line 
11 - 62 is of constant radius and is therefore one segment. Again, 
the coordinates of these nodes have been specified previously. iene 
62 - 113 is of constant angle, meeting the criterion for a segment. 
The coordinates of node 113 must be specified. Line 113 - 14/7 is 


of neither constant angle nor constant radius. ituaerebore., Le cannot 
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be treated as a segment and the coordinates of all nodes along ise 
line must be sneluded in the Input data. In order to further control 
the departure angle of the Lines mentioned above, line 60 - 145 

7s) also included as a segment. The coordinates of the end nodes 

of this segment have been computed in the problem already and do 

not need to be jneludcdmen che 1npue Gata. As with Line Sele 
the exact shape of line 60 - 145 is not important. Similarly, on 
the back face, lines SGi= Solon 826 = 941, and 864 - 949 are treated 
as segments. Although the nodes along line 81> - 951 are at the 
origin, care must be taken to insure that the.coordinates of these 
nodes are specified at appropriate ang les in such a way that this 
point will have the same "shape" as the front face. The radial 
component of anese coordinates is zero. This care is necessary 
because all computations in the problem are made with ew Tei repe seat | 
coordinates. Were the coordinates of the nodes along this Mab ete: 
specified at a radius of zero and an angle of zero, the mesh would 
be distorted. The following lines can also be treated as segments: 
815 - 866 and SGGue Oli oa ae coordinates of nodes 0255) eos and 

940 must be included in the input data for the same reason 46 was 
given for nodes 119, 139, and 136. 

To complete the process of specifying the coordinates of nodes 
alone the structure edges, the following lines will be treated as 
segments in the (‘-direction: I Oye eve 856 137 - 941, 60 - 864, 
w= 881, 94 - 898, 11 - 915, 11 - 815, 62 - 866, 113 - 917, and 
130 - 934. 

The procedure for determining the values of NL, Ml, and Ll 


for the segment division data willenee be discussed here. ie ats 
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felt that the explanation in Chapter 4, combined with the listing 
of the input cards for this problem, is Sufficient ly eebeam. 
The card containing the element field data is prepared after 
the above steps have been accomplished. Loading information and the 
request for punched cards should not be included in the input data 
until it has been determined that the mesh will be satisfactory. 
Approximately 45 seconds of computer time were required to solve 
this problem with QUAMEG. That time includes the 31 seconds required 
to compile the program. Figures UZpalen. 14a, and 16 are the plots 
of the structure mesh that was generated by QUAMEG. Table I is the 
listing of the required input cards described above. Figure l/ is 
the result of a 5 x 5 x 3 mesh generated with CUMEG for the same 


structure. Intermediate slices and the back surface are not shown. 
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Generated Mesh For Front Face 
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Figure 15. 





Generated Mesh For Third Intermediate Slice 
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APPENDIX B 
A CLASSICAL PROBLEM SOLVED WITH TRISOP 

TRISOP is a relatively new computer system which uses isoparametric 
finite elements to obtain stresses, strains, and displacements of 
three-dimensional structures that are subjected to static loads. 

Although a number of test problems have been conducted with 
TRISOP, the input preparation was too voluminous to contemplate 
a full evaluation of the system. Without QUAMEG, the difficulties 
discussed in Chapter 1 limited the size and number of problems that 
could be solved. In order to evaluate TRISOP as a complete system, 
a classical problem was selected and a convergence study was made, 
using QUAMEG and TRISOP. 

The problem selected was that of the pinched, thin-shell cylinder. 
Obviously, this problem should be solved using thin-shell elements 
for maximum efficiency. The purpose of this test was to determine 
if the three-dimensional solution would converge to the same results 
as are reported in Ref. 3 and also to determine how well thin-shell 
action is captured by TRISOP. The results Jenene by Cantin [3] 
are included in Table II with the results obtained in-the test. The 


problem is described in Figure 18. 


Conclusions: 
Although the results of this convergence study indicate a 
tendency to approach the same results obtained by thin-shell theory, 


acceptable results can be obtained oniy with a very fine mesh. 
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TRISOP can be used to solve problems for which there are no 
closed-form or accurate numerical methods of solution available. 
However, the user must be prepared to pay the price of long compu- 
tational times, for some problems. 


Table II 
Displacement at the applied load for a pinched cylinder. 


Bogner et al. Cantin 
No. No 
of Mesh Displace., of Mesh Displac., 
ale ens Gn 19h — 
72 he Som 2 -0.0802 36 ero -0,0921 
96 1x 3 -0.1026 48 aes 3 -0.1072 
cZ0 lee -0. 1087 60 ete ei -0.1099 
108 2x 2 -0.0808 54 Pome? -0.0931 
144 2) xs =e OCS 72 2x 3 -0.1085 
180 Ze -0.1098 90 aoe -0.1113 
150 4x4 -0.1126 
294 6x 6 -0.1137 
486 8 x 8 -0.1139 
726 10 x 10 -0.1139 
Hanson 
No. : 
of Mesh Displac., 
Eq. | al ge 
le 2 x 2 -0.00177 
243 Ji, SANE NP -0.00178 
423 EE IE: -0.00196 
465 ee Se cy -0.00375 
945 lex 6 x26 -0.03002 
1593 oxo xo -0.08702 
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APPENDIX C 
COMPUTER LISTING (QUAMEG) 
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COMPUTE ELEMENT CONNECTIVITY 
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